perm filename SAITR2.FAI[SYS,BGB] blob
sn#001408 filedate 1972-01-08 generic text, type T, neo UTF8
00100
00200 COMMENT ⊗
00300
00400 THIS IS THE FILE OF SAIL NUMERIC ROUTINES.
00500 IT MAY BE COMPILED STAND-ALONE, IN WHICH CASE
00600 ALL THE ENTRIES ARE MERELY DECLARED 'INTERNAL'.
00700
00800 ALTERNATIVELY, IT MAY BE USED TO PUT TOGETHER
00900 A LIBRARY. THE VARIABLE TELLME IS USED TO DETERMINE
01000 WHICH LIBRARY IS BEING MADE!!
01100
01200 ⊗
01300
01400
01500
01600
01700 TITLE NUMERC
01800 SUBTTL NUMERICAL ROUTINES FOR SAIL
01900
02000
02100 ↓P←17
02200
00100 BEGIN UNDER
00200 INTERNAL UNDERFLOW
00300 EXTERNAL JOBTPC,JOBAPR
00400 OPDEF ERROR [5B8]
00500
00600 ;The parameter to underflow controls the various messages:
00700 ;the 01 bit turns on floating underflow,
00800 ;the 02 bit turns on floating overflow
00900 ;the 04 bit truns on zero divide
01000 ;the 10 bit turns on fixed overflow
01100
01200 UNDERF: MOVEI 2,10
01300 JFCL 17,.+1; CLEAR ANY PREVIOUSLY SET FLAGS
01400 SKIPN 1,-1(17); CHANGED FROM SKIPE BECAUSE SENSE OF FLAG WAS WRONG
01500 SETZ 2,
01600 MOVE 3,[XWD -4,UNDFLG]
01700 LP: SETZM (3)
01800 TRNE 1,1
01900 SETOM (3) ;ENABLE ERROR MESSAGE ON BIT=1
02000 LSH 1,-1
02100 AOBJN 3,LP
02200 SUB 17,[(2)2]
02300 MOVEI 1,FLTOV
02400 MOVEM 1,JOBAPR
02500 CALLI 2,16 ;SET APR FLAGS
02600 MOVE 1,2(17)
02700 TLZ 1,440140 ;CLEAR PREV FLAGS
02800 JRST 2,@1 ;RETURN
02900
03000 FLTOV: MOVEM 1,SAVE1
03100 MOVE 1,JOBTPC
03200 MOVEM 1,OVPCWD
03300 TLNN 1,100 ; is it a floating underflow?
03400 JRST OV ;no
03500 MOVE 1,-1(1) ;get opcode which caused it
03600 TLNN 1,40000 ;test for standard flt pt opcode
03700 TLZ 1,2000 ;change for FSC
03800 DPB 1,[POINT 29,.+2,35] ;modify the SETZ
03900 MOVE 1,SAVE1 ;restore ACs
04000 SETZ 0, ;zero ac and/or memory
04100 MOVEM 1,SAVE1
04200 JSR JFCLTST
04300 SKIPN UNDFLG
04400 JRST WO ;dont print message
04500 MOVE 1,BP1
04600 JSR NUMOUT
04700 OUTSTR MESS1
04800 WO: MOVE 1,JOBTPC
04900 TLZ 1,440140 ;zero the error bits
05000 MOVEM 1,JOBTPC
05100 MOVE 1,SAVE1
05200 JRST 2,@JOBTPC ;return
05300
05400 OV: JSR JFCLTST
05500 TLNN 1,40000 ;was it a floating overflow?
05600 JRST ZDIV ;no
05700 SKIPN FOVFLG
05800 JRST WO ;dont print flt over message
05900 MOVE 1,BP2
06000 JSR NUMOUT
06100 ERROR 1,MESS2
06200 JRST WO
06300
06400 ZDIV: TLNN 1,40 ;zero divide?
06500 JRST NOTIN ;no
06600 SKIPN ZDFLG
06700 JRST WO ;dont print zero divide message
06800 MOVE 1,BP4
06900 JSR NUMOUT
07000 ERROR 1,MESS4
07100 JRST WO
07200
07300 NOTIN: SKIPN OVFLG
07400 JRST WO ;dont print fixed pt overflow message
07500 MOVE 1,BP3
07600 JSR NUMOUT
07700 ERROR 1,MESS3
07800 JRST WO
07900
08000 JFCLTST:0
08100 MOVE 1,JOBTPC
08200 MOVEM 2,SAVE2
08300 HLRZ 2,(1) ; get next opcode
08400 CAIE 2,(<JFCL>)
08500 JRST [MOVE 2,SAVE2 ; not JFCL, return
08600 JRST @JFCLTST]
08700 MOVE 2,SAVE2 ; a JFCL
08800 TRNN 1,-1
08900 JRST [MOVE 1,SAVE1 ;address of JFCL is zero, execute the JFCL
09000 JRST @JOBTPC]
09100 HRRM 1,JLOC ;go to address in JFCL
09200 MOVE 1,SAVE1
09300 JLOC: JRST @
09400
09500 NUMOUT: 0
09600 MOVEM 1,PTR
09700 MOVEM 2,SAVE2
09800 MOVEI 2,6
09900 MOVE 1,JOBTPC
10000 HRLZI 1,-1(1)
10100 L1: ROT 1,3
10200 IORI 1,60
10300 IDPB 1,PTR
10400 HLRI 1,
10500 SOJG 2,L1
10600 MOVE 2,SAVE2
10700 JRST @NUMOUT
10800
10900 INTERNAL OVPCWD
11000 ↑OVPCWD: 0
11100 PTR: 0
11200 UNDFLG: 0
11300 FOVFLG: 0
11400 ZDFLG: 0
11500 OVFLG: 0
11600 SAVE1: 0
11700 SAVE2: 0
11800 BP1: POINT 7,MESS1+6,20
11900 BP2: POINT 7,MESS2+6,13
12000 BP3: POINT 7,MESS3+4,20
12100 BP4: POINT 7,MESS4+5,13
12200 MESS1: ASCIZ /FLOATING UNDERFLOW OCCURED, PC = 000000
12300 /
12400 MESS2: ASCIZ/FLOATING OVERFLOW OCCURED, PC = 000000/
12500 MESS3: ASCIZ/OVERFLOW OCCURED, PC = 000000/
12600 MESS4: ASCIZ /ZERO DIVIDE OCCURED, PC = 000000/
12700
12800 BEND UNDER
00100 BEGIN $EXP3
00200 ;FROM V.021 22-SEP-69
00300 ;----------------------------------------------------
00400 ;SINGLE PRECISION FORTRAN IV EXP.3 FUNCTIONS
00500 ;THESE ROUTINES CALCULATE A FLOATING POINT NUMBER RAISED TO A
00600 ;FLOATING POINT POWER. THE CALCULATION IS
00700 ; A**B= EXP(B*LOG(A))
00800
00900 ;IF THE EXPONENT IS AN INTEGER LESS THAN 2**35 IN MAGNITUDE, THE
01000 ;RESULT WILL BE COMPUTED USING "$EXP2" AND THE ANSWER
01100 ;WILL HAVE THE CORRECT SIGN. (REMEMBER THAT THE "INTEGER"
01200 ;HAS ONLY 27 SIGNIFCANT BITS.)
01300 ;SINCE NEGATIVE NUMBERS RAISED TO NON-INTEGER POWERS YIELD
01400 ;COMPLEX ANSWERS, THE MAIN ALGORITHM CALCULATES
01500 ; EXP(B*LOG(ABSF(A))) = REALPART(A↑B).
01600
01700 ;THUS, SPECIAL CASES ARE:
01800 ; A↑0 = 1.
01900 ; 0↑B, B GEQ 0, = 0.
02000 ; 0↑B, B LESS THAN 0, = INF.
02100 ; A↑B, A LESS THAN 0,B INTEGER, = (-1)↑B*ABS(A)↑B, USES $EXP2.
02200 ; A↑B, A LESS THAN 0,B NON-INTEGRAL, = REALPART(A↑B).
02300
02400
02500 ;THE CALLING SEQUENCE FOR THE ROUTINE IS AS FOLLOWS:
02600 ; PUSH P,BASE
02700 ; PUSH P,EXPNT
02800 ; PUSHJ P,$EXP3
02900 ;THE RESULT IS RETURNED IN AC 1.
03000
03100
03200 ;ACCUMULATOR DEFINITIONS
03300 A←15
03400 B←1
03500 C←13
03600 D←14
03700
03800 INTERNAL EXP3,EXP3.0
03900
04000 EXP3.0: PUSH P,15
04100 PUSH P,0
04200 PUSH P,1
04300 PUSHJ P,EXP3
04400 POP P,15
04500 MOVE 0,1
04600 POPJ P,
04700
04800 EXP3:
04900 $EXP3: MOVE A,-2(P) ;AC A ← BASE
05000 MOVE B,-1(P) ;AC B ← EXPONENT
05100 JUMPE B,[MOVSI B,(1.0);BASE**0, RETURNS 1
05200 JRST EXIT]
05300 JUMPN A,EXP30A ;GO AHEAD IF BASE NE 0.
05400 JUMPGE B,SET0 ;EXIT IF BASE = 0, EXP GREATER THAN= 0,
05500 OUTSTR [ASCIZ /EXP3: 0 TO A NEGATIVE POWER; LARGEST LEGAL NUMBER RETURNED/]
05600 HRLOI B,377777 ;ANS.=+INFINITY
05700 JRST EXIT ;AND EXIT.
05800 SET0: SETZ B,0 ;ANSWER = 0.
05900 JRST EXIT ;RETURN.
06000
06100 EXP30A:
06200 MOVM D,B ;SET EXP. POSITIVE.
06300 MOVEI C,0 ;CLEAR AC C TO ZERO
06400 LSHC C,11 ;SHIFT 9 PLACES LEFT
06500 SUBI C,200 ;TO OBTAIN SHIFTING FACTOR
06600 JUMPLE C,EXP3GO ;IS C GREATER THAN 0 ,IF N0, NOT AN INTEGER.
06700 HRR B,C ;SET UP B AS AN INDEX REG.
06800 MOVEI C,0 ;CLEAR OUT AC C
06900 LSH D,-1 ;RIGHT ADJUST EXP TO BIT 1.
07000 ASHC C,(B) ;SHIFT LFT BY CONTENTS OF B.
07100 ;THUS PUT INTEGER PART IN C,
07200 ;SET OVERFLOW IF ANY BITS FALL
07300 ;OFF THE LEFT END OF C.
07400 JFCL EXP3GO ;IF OVERFLOW, GO TO EXP3GO.
07500 ;THE SPECIAL INTERRUPT CODE
07600 ;WILL CAUSE THIS BRANCH TO BE
07700 ;TAKEN IF OVERFLOW OCCURS.
07800 ;NO RESULTS WILL BE CHANGED.
07900 ;EXPONENT IS INTEGRAL, BUT
08000 ;TOO LARGE TO TREAT AS SUCH.
08100 ;(D MUST BE EMPTY.)
08200 JUMPN D,EXP3GO ;IS EXPONENT AN INTEGER ?
08300 SKIPGE -1(P) ;YES, WAS IT NEG. ?
08400 MOVNS C ;YES, NEGATE IT
08500 PUSH P,A
08600 PUSH P,C
08700 PUSHJ P,$EXP2 ;OBTAIN RESULT USING $EXP2
08800 JRST EXIT ;RETURN
08900 EXP3GO:
09000 JUMPGE A,GETLOG ;GO TO TAKE LOG OF A.
09100 OUTSTR [ASCIZ /EXP3:NEGATIVE QUANTITY TO A NON-INTEGRAL POWER;
09200 REALPART OF RESULT RETURNED./];
09300
09400 MOVMM A,A ;GET ABS(BASE) (BASE NE 0 OR 1.
09500 GETLOG: PUSH P,A ;CALCULATE LOG OF A
09600 PUSHJ P,$LOG
09700 FMPR B, -1(P) ;CALCULATE B*LOG(A)
09800 PUSH P,B ; CALCULATE EXP(B*LOG(A))
09900 PUSHJ P,$EXP
10000 EXIT:
10100 SUB P,[XWD 3,3]
10200 JRST @3(P)
10300
10400
10500 BEND $EXP3
00100 ;-------------------------------------------------------
00200 ;SINGLE PRECISION EXP.2 FUNCTIONS
00300 ;THESE ROUTINES CALCULATE A FLOATING POINT NUMBER TO A FIXED
00400 ;POINT POWER. THE CALCULATION IS A**B, WHERE B IS OF THE FORM
00500
00600 ; B=Q(0) + Q(1)*2 + Q(2)*4 + ...WHERE Q(I)=0 OR 1
00700
00800 ;THERE ARE NO RESTRICTIONS ON THE BASE OR EXPONENT
00900
01000 ;THE CALLING SEQUENCE FOR THE ROUTINE IS THE FOLLOWING:
01100 ; PUSH P,BASE
01200 ; PUSH P,EXPNT
01300 ; PUSHJ P,$EXP2
01400 ;THE RESULT IS RETURNED IN AC 1.
01500
01600
01700 BEGIN $EXP2
01800 EXTERN OVPCWD
01900
02000 ;ACCUMULATOR DEFINITIONS
02100 B←1
02200 A←13
02300 C←14
02400 D←15
02500
02600 INTERNAL EXP2,EXP2.0
02700
02800 EXP2.0: PUSH P,15
02900 PUSH P,0
03000 PUSH P,1
03100 PUSHJ P,EXP2
03200 POP P,15
03300 MOVE 0,1
03400 POPJ P,
03500
03600 EXP2:
03700 ↑$EXP2: MOVE A,-2(P) ;AC A ← BASE.
03800 MOVE B,-1(P) ;AC B ← EXPONENT.
03900 JUMPE B,[MOVSI B,(1.0);BASE**0, RETURNS 1
04000 JRST EXIT]
04100 JUMPN A,EXP2A ;GO AHEAD IF BASE NE 0.
04200 JUMPGE B,SET0 ;EXIT IF BASE =0, EXP GREATER THAN 0,
04300 OUTSTR [ASCIZ /EXP2: 0 TO A NEGATIVE POWER; LARGEST LEGAL NUMBER RETURNED/];
04400 HRLOI B,377777 ;AN ANSWER OF INFINITY.
04500 JRST EXIT ;RETURN
04600
04700 EXP2A:
04800 MOVE D, B ;PUT EXPONENT IN D.
04900 MOVSI B, (1.0) ;GET 1.0 IN ACCUMULATOR B.
05000 MOVM C, D ;MAKE EXPONENT POSITIVE IN C.
05100 JFCL MININF ;IF EXP WAS 400000,,0 GO TO MININF.
05200 JRST FEXP2 ;ENTER MAIN PART OF PROGRAM.
05300
05400 FEXP1: FMPR A, A ;FORM A**N, FLOATING POINT.
05500 JFCL OVER ;IF OVER/UNDERFLOW, GO TO OVER.
05600 LSH C, -1 ;SHIFT EXPONENT FOR NEXT BIT.
05700 FEXP2: TRZE C, 1 ;IS THE BIT ON?
05800 FMPR B, A ;YES, MULTIPLY ANSWER BY A**N.
05900 JFCL OVER ;IF OVER/UNDERFLOW, GO TO OVER.
06000 JUMPN C, FEXP1 ;UPDATE A**N UNLESS ALL THROUGH.
06100 INVTST: JUMPGE D, EXIT ;JUMP IF RECIPROCAL IS NOT NEEDED.
06200 INV: MOVSI C, (1.0) ;GET 1.0 IN C.
06300 FDVRM C, B ;FORM 1/(A**B) FOR NEG. EXPONENT.
06400 JFCL MSGINF ;IF 1/0, SPECIAL INTERRUPT CODE
06500 ;WILL SET B TO +INF, AND WE'LL BE
06600 ;SENT TO MSGINF. (1/+-INF IS OK
06700 ;AND WILL NOT CAUSE AN INTERRUPT.)
06800
06900 EXIT:
07000 SUB P,[XWD 3,3]
07100 JRST @3(P) ;RETURN.
07200
07300
07400 OVER:
07500 MOVE C,OVPCWD ;PICK UP FLAGS.(SET BY INTERRUPT CODE)
07600 TLNN C,(1B11) ;SKIP IF UNDERFLOW.
07700 JRST OFM ;JUMP IF OVERFLOW.
07800 JUMPL D,SETBIG ;JUMP IF EXP LESS THAN 0.
07900 MSG0: OUTSTR [ASCIZ /EXP2:RESULT MAGNITUDE TOO SMALL;0 RETURNED./];
08000 SET0: SETZ B,0 ;SET B TO 0.
08100 JRST EXIT ;RETURN
08200
08300 SETBIG: HRLOI B,377777 ;RESULT IS +- 1/0, WITH
08400 ;SIGN=SIGN(BASE)↑PARITY(EXP).
08500 SKIPL -2(P) ;IS A NEGATIVE?
08600 JRST MSGINF ;NO
08700 TRNE D,1 ;ANS. IS LESS THAN 0, IF A LESS THAN 0 AND
08800 MOVNS B ;THE EXP. WAS ODD.
08900 MSGINF: OUTSTR [ASCIZ /EXP2:RESULT MAGNITUDE TOO LARGE; BIG SIGNED # RETURNED./];
09000 JRST EXIT ;RETURN
09100
09200 OFM: JUMPG D,SETBIG ;IF EXP GREATER THAN 0,RESULT IS +-INF.
09300 JRST MSG0 ;RESULT IS 1/+-INF =0.
09400
09500 MININF: HRLOI C,377777 ;SET EXP = +INFINITY.
09600 MOVE B,A ;SET B UP BY 1 FACTOR OF A.
09700 JRST FEXP2 ;GO TO MAIN ROUTINE.
09800
09900
10000 SAVEA: 0 ;TEMP FOR A.
10100 SAVEB: 0 ;TEMP FOR B.
10200 SAVEC: 0 ;TEMP FOR C.
10300
10400 BEND $EXP2
00100 ;ORIGINALLY: ALOG V.027 FROM V.022 18-DEC-69 FROM V.020
00200 ; 17-JUL-70 /KK/DMN
00300 ;---------------------------------------------------------
00400 ;FLOATING POINT SINGLE PRECISION LOGARITHM FUNCTION
00500 ;LOG(ABSF(X)) IS CALCULATED BY THE SUBROUTINE, AND AN
00600 ;ARGUMENT OF ZERO IS RETURNED AS MINUS INFINITY.
00700
00800 ;$LOG IS THE ENTRY POINT FOR LOGE(X), AND
00900 ;$LOG10 IS THE ENTRY POINT FOR LOG10(X).
01000 ;FOR LOGE(X), THE ALGORITHM IS:
01100 ; LOGE(X) = (I + LOG2(F))*LOGE(2)
01200 ; WHERE X = (F/2)*2↑(I+1), AND LOG2(F) IS GIVEN BY
01300 ; LOG2(F) = C1*Z + C3*Z↑3 + C5*Z↑5 - 1/2
01400 ; AND Z = (F-SQRT(2))/(F+SQRT(2))
01500 ;FOR LOG10(X), THE ALGORITHM IS:
01600 ; LOG10(X) = LOGE(X)*LOG10(E)
01700
01800 ;THE CALLING SEQUENCE FOR THE ROUTINE IS AS FOLLOWS:
01900 ; PUSH 17,ARG
02000 ; PUSHJ 17,$LOG OR $LOG10
02100 ;THE ANSWER IS RETURNED IN ACCUMULATOR 1
02200
02300
02400 BEGIN $LOG
02500
02600 A←13
02700 B←14
02800 C←15
02900 D←1
03000
03100 INTERNAL LOG,LOG10
03200
03300 LOG10:
03400 $LOG10: ;ENTRY TO LOG TO THE BASE 10 ROUTINE.
03500 SKIPA D,LOG10A ;SETUP FOR LOG10 FACTOR.
03600
03700 LOG:
03800 ↑$LOG: ;ENTRY TO LOG TO THE BASE E ROUTINE.
03900 MOVSI D, (1.0) ;SET D TO 1.0.
04000 SKIPGE 0, -1(P) ;SKIP IF ARG NOT NEGATIVE.
04100 OUTSTR [ASCIZ /LOG:NEGATIVE ARGUMENT. REAL PART RETURNED./]
04200
04300 MOVM A,-1(P) ;GET ABSF(X)
04400 JUMPE A, LZERO ;CHECK FOR ZERO ARGUMENT
04500 CAMN A, D ;CHECK FOR 1.0 ARGUMENT
04600 JRST ZERANS ;IT IS 1.0 RETURN ZERO ANS.
04700 ASHC A, -33 ;SEPARATE FRACTION FROM EXPONENT
04800 ADDI A, 211000 ;FLOAT THE EXPONENT AND MULT. BY 2
04900 MOVSM A, C ;NUMBER NOW IN CORRECT FL. FORMAT
05000 MOVSI A, 567377 ;SET UP -401.0 IN A
05100 FADM A, C ;SUBTRACT 401 FROM EXP.*2
05200 ASH B, -10 ;SHIFT FRACTION FOR FLOATING
05300 TLC B, 200000 ;FLOAT THE FRACTION PART
05400 FAD B, L1 ;B = B-SQRT(2.0)/2.0
05500 MOVE A, B ;PUT RESULTS IN A
05600 FAD A, L2 ;A = A+SQRT(2.0)
05700 FDV B, A ;B = B/A
05800 MOVEM B, LZ ;STORE NEW VARIABLE IN LZ
05900 FMP B, B ;CALCULATE Z↑2
06000 MOVE A, L3 ;PICK UP FIRST CONSTANT
06100 FMP A, B ;MULTIPLY BY Z↑2
06200 FAD A, L4 ;ADD IN NEXT CONSTANT
06300 FMP A, B ;MULTIPLY BY Z↑2
06400 FAD A, L5 ;ADD IN NEXT CONSTANT
06500 FMP A, LZ ;MULTIPLY BY Z
06600 FAD A, C ;ADD IN EXPONENT TO FORM LOG2(X)
06700 FMP A, L7 ;MULTIPLY TO FORM LOGE(X)
06800 FMPR D, A ;GET LOGE OR LOG10 IN D.
06900 EXIT: SUB P,[XWD 2,2]
07000 JRST @2(P)
07100 LZERO: OUTSTR [ASCIZ /LOG: ARGUMENT = 0; MINUS INFINITY RETURNED/];
07200 MOVE D,MIFI ;PICK UP MINUS INFINITY
07300 JRST EXIT ;EXIT
07400 ZERANS: MOVEI D, 0 ;MAKE ANSWER ZERO
07500 JRST EXIT
07600
07700 ;CONSTANTS
07800
07900 LOG10A: 177674557305
08000 L1: 577225754146 ;-0.707106781187
08100 L2: 201552023632 ;1.414213562374
08200 L3: 200462532521 ;0.5989786496
08300 L4: 200754213604 ;0.9614706323
08400 L5: 202561251002 ;2.8853912903
08500 L7: 200542710300 ;0.69314718056
08600 MIFI: 400000000001 ;LARGEST NEGATIVE FLOATING NUMBER
08700
08800 LZ: 0
08900
09000 BEND $LOG
00100 ;---------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION EXPONENTIAL FUNCTION
00300 ;IF X LESS THAN=-89.415..., THE PROGRAM RETURNS ZERO AS THE ANSWER
00400 ;IF X GREATER THAN=88.029..., THE PROGRAM RETURNS 377777777777 AS THE ANSWER
00500 ;THE RANGE OF THE ARGUMENT IS REDUCED AS FOLLOWS:
00600 ;EXP(X) = 2**(X*LOG(E)BASE2) = 2**(M+F)
00700 ;WHERE M IS AN INTEGER AND F IS A FRACTION
00800 ;2**M IS CALCULATED BY ALGEBRAICALLY ADDING M TO THE EXPONENT
00900 ;OF THE RESULT OF 2**F. 2**F IS CALCULATED AS
01000
01100 ;2**F = 2(0.5+F(A+B*F↑2 - F-C(F↑2 + D)**-1)**-1
01200
01300 ;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
01400 ; PUSH 17,ARG
01500 ; PUSHJ 17,$EXP
01600 ;THE ANSWER IS RETURNED IN ACCUMULATOR 1
01700
01800 BEGIN $EXP
01900 A←1
02000 B←13
02100 C←14
02200 D←15
02300
02400
02500 INTERNAL EXP
02600
02700
02800 EXP:
02900 ↑$EXP: ;ENTRY TO EXPONENTIAL ROUTINE
03000 MOVE B,-1(P) ;PICK UP THE ARGUMENT IN B
03100 CAMGE B,E77 ;IS EXP. LESS THAN -89.41...?
03200 JRST OUT2 ;YES, GO TO EXIT.
03300 CAMG B,E7 ;IS EXP. GREATER THAN +88.029...?
03400 JRST EXP1 ;GO TO STANDARD ALGORITHM.
03500 OUTSTR [ASCIZ /EXP: RESULT TOO LARGE; LARGEST LEGAL NUMBER RETURNED/];
03600 HRLOI 1, 377777 ;GET LARGEST FLOATING NUMBER
03700 JRST EXIT ;EXIT
03800
03900 OUT2: OUTSTR [ASCIZ /EXP: RESULT TOO SMALL; 0 RETURNED/];
04000 MOVEI 1,0 ;ANSWER IS 0.
04100 JRST EXIT ;EXIT
04200
04300
04400 EXP1:
04500 SETZM ES2 ;INITIALIZE ES2
04600 MULI B, 400 ;SEPARATE FRACTION AND EXPONENT
04700 TSC B, B ;GET A POSITIVE EXPONENT
04800 MUL C, E5 ;FIXED POINT MULTIPLY BY LOG2(E)
04900 ASHC C, -242(B) ;SEPARATE FRACTION AND INTEGER
05000 AOSG C ;ALGORITHM CALLS FOR MULT. BY 2
05100 AOS C ;ADJUST IF FRACTION WAS NEGATIVE
05200 HRRM C, EX1 ;SAVE FOR FUTURE SCALING
05300 JUMPG D,ASHH ;GO AHEAD IF ARG GREATER THAN 0.
05400 TRNN D,377 ;ARE ALL THESE BITS 0?
05500 JRST ASHH ;YES, GO AHEAD.
05600 ADDI D,200 ;NO, FIX UP.
05700 ASHH: ASH D, -10 ;MAKE ROOM FOR EXPONENT
05800 TLC D, 200000 ;PUT 200 IN EXPONENT BITS
05900 FADB D, ES2 ;NORMALIZE, RESULTS TO D AND ES2
06000 FMP D, D ;FORM X↑2
06100 MOVE A, E2 ;GET FIRST CONSTANT
06200 FMP A, D ;E2*X↑2 IN A
06300 FAD D, E4 ;ADD E4 TO RESULTS IN D
06400 MOVE B, E3 ;PICK UP E3
06500 FDV B, D ;CALCULATE E3/(F↑2 + E4)
06600 FSB A, B ;E2*F↑2-E3(F↑2 + E4)**-1
06700 MOVE C, ES2 ;GET F AGAIN
06800 FSB A, C ;SUBTRACT FROM PARTIAL SUM
06900 FAD A, E1 ;ADD IN E1
07000 FDVM C, A ;DIVIDE BY F
07100 FAD A, E6 ;ADD 0.5
07200 EX1: FSC A, 0 ;SCALE THE RESULTS
07300 EXIT: SUB P,[XWD 2,2]
07400 JRST @2(P)
07500
07600 E1: 204476430062 ;9.95459578
07700 E2: 174433723400 ;0.03465735903
07800 E3: 212464770715 ;617.97226953
07900 E4: 207535527022 ;87.417497202
08000 E5: 270524354513 ;LOG(E), BASE 2
08100 E6: 0.5
08200 E7: 207540074636 ;88.029...
08300 E77: 570232254037 ;-89.415986
08400 ES2: 0
08500
08600 BEND $EXP
00100 ;--------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION SINE AND COSINE FUNCTION
00300
00400 ;IF THE ARGUMENT IS IN DEGREES, THE PROPER ENTRY POINTS ARE
00500 ;$SIND AND $COSD, WHILE IF THE ARGUMENT IS IN RADIANS, THE
00600 ;PROPER ENTRY POINTS ARE $SIN AND $COS.
00700 ;$COSD CALLS $SIND TO CALCULATE SIND(PI/2+X)
00800 ;$COS CALLS $SIN TO CALCULATE SIN (PI/2+X)
00900 ;$SIND CALLS $SIN AFTER A CONVERSION FROM DEGREES TO RADIANS.
01000
01100 ;THIS ROUTINE CALCULATES SINES AFTER REDUCING THE ARGUMENT TO
01200 ;THE FIRST QUADRANT AND CHECKING THE OVERFLOW BITS TO DETERMINE
01300 ;THE QUADRANT OF THE ORIGINAL ARGUMENT.
01400 ;000 - 1ST QUADRANT
01500 ;001 - 2ND QUADRANT
01600 ;010 - 3RD QUADRANT
01700 ;011 - 4TH QUADRANT
01800 ;THE ALGORITHM USES A MODIFIED TAYLOR SERIES TO CALCULATE
01900 ;THE SINE OF THE NORMALIZED ARGUMENT.
02000
02100 ;THE ROUTINES ARE CALLED IN THE FOLLOWING MANNER:
02200 ; PUSH P,ARG
02300 ; PUSHJ P,$SIN (OR $COS,$SIND, OR $COSD)
02400 ;THE ANSWER IS RETURNED IN ACCUMULATOR 1
02500 BEGIN $SIN
02600
02700 A←13
02800 B←14
02900 C←15
03000
03100 INTERNAL SIN,SIND,COS,COSD
03200
03300 COSD:
03400 $COSD: ;ENTRY TO COSINE DEGREES ROUTINE.
03500 MOVE B,-1(P) ;PICK UP THE ARG.
03600 FADR B,CD1 ;ADD 90 DEGREES.
03700 JRST CONVER ;CONVERT TO RADIANS
03800 ;THEN GO TO SIN ROUTINE
03900
04000 SIND:
04100 $SIND: ;ENTRY TO SINE DEGREES ROUTINE.
04200 MOVE B,-1(P) ;PICK UP THE ARG.
04300 CONVER: FDVR B,SCD1 ;CONVERT TO RADIANS
04400 JFCL ;SUPPRESS ERROR MESSAGE ON UNDERFLOW.
04500 ;SPECIAL INTERRUPT CODE WILL SET
04600 ; B TO 0 ON UNDERFLOW
04700 JRST S1 ;ENTER SINE ROUTINE.
04800
04900 COS:
05000 $COS: ;ENTRY TO COSINE RADIANS ROUTINE.
05100 MOVE B,-1(P) ;PICK UP THE ARG.
05200 FADR B,PIOT ;ADD PI/2.
05300 JRST S1 ;ENTER SINE ROUTINE.
05400
05500
05600 SIN:
05700 $SIN: ;ENTRY TO SINE RADIANS ROUTINE.
05800 MOVE B,-1(P) ;PICK UP THE ARG.
05900 S1: MOVEM B,-1(P) ;SAVE THE ARG.
06000 MOVMS B ;GET ABS OF ARG.
06100 CAMG B,SP2 ;SIN(X)=X IF X LESS THAN 2↑-9.
06200 JRST S3A ;EXIT WITH ARG. IN B.
06300 FDV B,PIOT ;DIVIDE X BY PI/2.
06400 CAMG B,ONE ;IS X/(PI/2) LESS THAN 1.0 ?
06500 JRST S2 ;YES,ARG IN 1ST QUADRANT ALREADY.
06600 MULI B,400 ;NO,SEPARATE FRACTION AND EXP.
06700 ASH C,-202(B) ;GET X MODULO 2PI.
06800 JFCL ;SUPRESS ERROR MESSAGE FROM OVTRAP.
06900 ;SPECIAL INTERRUPT CODE WILL
07000 ;RETURN WITHOUT ATTEMPTING A
07100 ;FIXUP
07200 MOVEI B,200 ;PREPARE FLOATING FRACTION.
07300 ROT C,3 ;SAVE THREE BITS TO DETERMINE QUADRANT.
07400 LSHC B,33 ;ARGUMENT NOW IN THE RANGE (-1,1).
07500 FAD B,SP3 ;NORMALIZE THE ARGUMENT.
07600 JUMPE C,S2 ;REDUCED TO 1ST QUAD IF BITS 000.
07700 TLCE C,1000 ;SUBTRACT 1.0 FROM ARG IF BITS ARE
07800 FSB B,ONE ;001 OR 011.
07900 TLCE C,3000 ;CHECK FOR FIRST QUADRANT, 001.
08000 TLNN C,3000 ;CHECK FOR THIRD QUADRANT, 010.
08100 MOVNS B ;001,010.
08200 S2: SKIPGE -1(P) ;CHECK SIGN OF ORIGINAL ARG.
08300 MOVNS B ;SIN(-X)=-SIN(X).
08400 MOVEM B,-1(P) ;STORE REDUCED ARG.
08500 FMPR B,B ;CALCULATE X↑X
08600 MOVE A,SC9 ;GET 1ST CONSTANT.
08700 FMP A,B ;MULTIPLY BY X↑2
08800 FAD A,SC7 ;ADD IN NEXT CONSTANT.
08900 FMP A,B ;MULTIPLY BY X↑2.
09000 FAD A,SC5 ;ADD IN NEXT CONSTANT.
09100 FMP A,B ;MULTIPLY BY X↑2.
09200 FAD A,SC3 ;ADD IN NEXT CONSTANT.
09300 FMP A,B ;MULTIPLY BY X↑2.
09400 FAD A,PIOT ;ADD IN LAST CONSTANT.
09500 S2B: FMPR A,-1(P) ;MULTIPLY BY X.
09600 MOVE 1,A ;ANSWER IN 1
09700 SKIPA 0,0 ;
09800 S3A: MOVE 1,-1(P) ;ANSWER IN 1.
09900 SUB P,[XWD 2,2]
10000 JRST @2(P)
10100
10200 SC3: 577265210372
10300 SC5: 175506321276
10400 SC7: 606315546346
10500 SC9: 164475536722
10600
10700 SP2: 170000000000
10800 SP3: 0
10900 CD1: 90.0
11000 SCD1: 206712273406
11100 PIOT: 201622077325
11200 ONE: 1.0
11300
11400
11500 BEND $SIN
00100 ;--------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION SQUARE ROOT FUNCTION
00300 ;THE SQUARE ROOT OF THE ABSOLUTE VALUE OF THE ARGUMENT IS
00400 ;CALCULATED. THE ARGUMENT IS WRITTEN IN THE FORM
00500 ; X= F*(2**2B) WHERE 0 LESS THAN F LESS THAN 1
00600 ;SQRT(X) IS THEN CALCULATED AS (SQRT(F))*(2**B)
00700 ;SQRT(F) IS CALCULATED BY A LINEAR APPROXIMATION, THE NATURE
00800 ;OF WHICH DEPENDS ON WHETHER 1/4 LESS THAN F LESS THAN 1/2 OR 1/2 LESS THAN F LESS THAN 1,
00900 ;FOLLOWED BY TWO ITERATIONS OF NEWTON'S METHOD.
01000
01100 ;THE CALLING SEQUENCE FOR THE SQUARE ROOT IS AS FOLLOWS:
01200 ; PUSH 17,ARG
01300 ; PUSHJ 17,$SQRT
01400 ;THE ANSWER IS RETURNED IN ACCUMULATOR 1.
01500
01600 BEGIN $SQRT
01700
01800 A←13
01900 B←14
02000
02100
02200 INTERNAL SQRT
02300
02400 SQRT:
02500 ↑$SQRT: ;ENTRY TO SQUARE ROOT ROUTINE
02600 SKIPG B,-1(P) ;PICK UP ARG. CHECK IF GREATER THAN 0
02700 JRST SQRT4 ;NO, HANDLE NON-POSITIVE ARGUMENT
02800
02900 MOVEI A,0 ;GET EXPONENT TO A
03000 LSHC A,=9
03100 SUBI A,201 ;GET TRUE EXPONENT + 1
03200 ROT A,-1 ;DIVIDE BY 2
03300 HRRM A,SQRT2 ;AND STORE FOR FLOATING SCALE INST.
03400 JUMPL A,SQRT3 ;JUMP IF FRACTION GREATER THAN .5
03500
03600 ;FRACTION LESS THAN .5
03700 LSH B,=-9 ;RESTORE POSITION OF FRACTION IN B
03800 FSC B,177 ;AND FIX UP EXPONENT .25 LESS THAN F LESS THAN .5
03900 MOVEM B,F ;SAVE FRACTION
04000 ;COMPUTE LINEAR APPROX #1
04100 FMPRI B,200640
04200 FADRI B,177465
04300
04400 SQRT1: MOVE A,F ;1ST ITERATION OF NEWTON
04500 FDV A,B ; F/APPROX
04600 FAD B,A ; APPROX + F/APPROX
04700 FSC B,-1 ; .5*( APPROX + F/APPROX)
04800 MOVE A,F ;2ND ITERATION OF NEWTON
04900 FDV A,B ; F/APPROX
05000 FADR A,B ; APPROX + F/APPROX
05100 SQRT2: FSC A,0 ;HALVE AND SCALE EXPONENT
05200 MOVE 1,A ;RESULT IN AC 1
05300 JRST EXIT ;RETURN
05400
05500 ;HERE ON F GREATER THAN= .5
05600 SQRT3: LSH B,=-9 ;RESTORE POSITION OF FRACTION IN B
05700 FSC B,200 ;AND FIX UP EXPONENT .5 LESS THAN= F LESS THAN 1
05800 MOVEM B,F ;SAVE FRACTION
05900 ;COMPUTE LINEAR APPROX #2
06000 FMPRI B,200450
06100 FADRI B,177660
06200 JRST SQRT1 ;NOW GO ITERATE
06300 SQRT4: JUMPE B,ZERO
06400 OUTSTR [ASCIZ /SQRT: NEGATIVE ARGUMENT - 0 RETURNED/];
06500 ZERO: MOVEI 1,0 ;HERE ON NON-POSITIVE ARG. RETURN ZERO
06600 EXIT: SUB P,[XWD 2,2]
06700 JRST @2(P)
06800
06900 F: BLOCK 1 ;STORE FRACTION HERE
07000
07100
07200 BEND $SQRT
00100 ;---------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION ARCTANGENT FUNCTION
00300 ;ATAN(X) = X(B0+A1(Z+B1-A2(Z+B2-A3(Z+B3)**-1)**-1)**-1)
00400 ;WHERE Z=X↑2, IF 0 LESS THAN X LESS THAN= 1
00500
00600 ;IF X GREATER THAN 1, THEN ATAN(X) = PI/2 - ATAN(1/X)
00700 ;AC D IS USED INTERNALLY TO KEEP TRACK OF CASES
00800 ;IF X GREATER THAN 1, THEN RH(D) =-1, AND LH(D) = -SGN(X)
00900 ;IF X LESS THAN 1, THEN RH(D) = 0, AND LH(D) = SGN(X)
01000
01100 ;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
01200 ; PUSH 17,ARG
01300 ; PUSHJ 17,$ATAN
01400 ;THE ANSWER IS RETURNED IN ACCUMULATOR 1
01500
01600 BEGIN $ATAN
01700
01800 A←1
01900 B←13
02000 C←14
02100 D←15
02200
02300 INTERNAL ATAN
02400
02500 ATAN:
02600 ↑$ATAN: ;ENTRY TO ARCTANGENT ROUTINE
02700 MOVE A,-1(P) ;PICK UP THE ARGUMENT IN A
02800 ATAN1: MOVM B, A ;GET ABSF OF ARGUMENT
02900 CAMG B, A1 ;IF X LESS THAN 2↑-33, THEN RETURN WITH...
03000 JRST EXIT ;ATAN(X) = X
03100 HLLO D, A ;SAVE SIGN, SET RH(D) = -1
03200 CAML B, A2 ;IF A GREATER THAN 2↑33, THEN RETURN WITH
03300 JRST AT4 ;ATAN(X) = PI/2
03400 MOVSI C, (1.0) ;FORM 1.0 IN C
03500 CAMG B, C ;IS ABSF(X) GREATER THAN 1.0?
03600 TRZA D, -1 ;IF B .LE. 1.0, THEN RH(D) = 0
03700 FDVM C, B ;B IS REPLACED BY 1.0/B
03800 TLC D, (D) ;XOR SIGN WITH .G. 1.0 INDICATOR
03900 MOVEM B, C3 ;SAVE THE ARGUMENT
04000 FMP B, B ;GET B↑2
04100 MOVE C, KB3 ;PICK UP A CONSTANT
04200 FAD C, B ;ADD B↑2
04300 MOVE A, KA3 ;ADD IN NEXT CONSTANT
04400 FDVM A, C ;FORM -A3/(B↑2 + B3)
04500 FAD C, B ;ADD B↑2 TO PARTIAL SUM
04600 FAD C, KB2 ;ADD B2 TO PARTIAL SUM
04700 MOVE A, KA2 ;PICK UP -A2
04800 FDVM A, C ;DIVIDE PARTIAL SUM BY -A2
04900 FAD C, B ;ADD B↑2 TO PARTIAL SUM
05000 FAD C, KB1 ;ADD B1 TO PARTIAL SUM
05100 MOVE A, KA1 ;PICK UP A1
05200 FDV A, C ;DIVIDE PARTIAL SUM BY A1
05300 FAD A, KB0 ;ADD B0
05400 FMP A, C3 ;MULTIPLY BY ORIGINAL ARGUMENT
05500 TRNE D, -1 ;CHECK .G. 1.0 INDICATOR
05600 FSB A, PIOT ;ATAN(A) = -(ATAN(1/A)-PI/2)
05700 SKIPA 0,0 ;SKIP
05800 AT4: MOVE A, PIOT ;GET PI/2 AS ANSWER
05900 SKIPGE D ;LH(D) = -SGN(B) IF B GREATER THAN 1.0
06000 MOVNS A ;NEGATE ANSWER
06100 EXIT: SUB P,[XWD 2,2]
06200 JRST @2(P)
06300 A1: 145000000000 ;2**-33
06400 A2: 233000000000 ;2**33
06500 KB0: 176545543401 ;0.1746554388
06600 KB1: 203660615617 ;6.762139240
06700 KB2: 202650373270 ;3.316335425
06800 KB3: 201562663021 ;1.448631538
06900 KA1: 202732621643 ;3.709256262
07000 KA2: 574071125540 ;-7.106760045
07100 KA3: 600360700773 ;-0.2647686202
07200 C3: 0
07300 PIOT: 201622077325 ;PI/2
07400
07500
07600 BEND $ATAN
00100 ;---------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION ARCSINE FUNCTION
00300 ;THE ARCSINE IS CALCULATED WITH THE FOLLOWING ALGORITHM:
00400
00500 ; ASIN(X) = ATAN(X/SQRT(1-X↑2))
00600
00700 ;THE RANGE OF DEFINITION FOR ASIN IS (-1.0,1.0)
00800 ;OTHER ARGUMENTS WILL CAUSE AN ERROR MESSAGE TO BE
00900 ;TYPED AND AN ANSWER OF ZERO TO BE RETURNED.
01000 ;CALLING SEQUENCE:
01100 ; PUSH P,ARG
01200 ; PUSHJ P,$ASIN
01300 ;
01400 ;RESULT RETURNED IN AC 1.
01500 ;
01600 ;
01700 BEGIN $ASIN
01800
01900 A←13
02000 B←1
02100
02200 INTERNAL ASIN
02300
02400 ASIN:
02500 $ASIN: ;ENTRY TO ASIN ROUTINE
02600 MOVM B,-1(P) ;GET MAGNITUDE OF ARG. IN B
02700 CAMLE B,ONE ;IS THE MAGNITUDE OF THE ARG. LE 1.0?
02800 JRST TOOLRG ;NO, GO TO ERROR RETURN.
02900 MOVN A,-1(P) ;GET THE NEGATIVE OF ARG
03000 FMP A,-1(P) ;CALCULATE -(X↑2)
03100 JFCL ;SUPPRESS ERROR MESSAGE FROM OVTRAP, IF NECESSARY.
03200 ;ON UNDERFLOW, THE SPECIAL INTERRUPT
03300 ;CODE SETS A TO 0
03400 FAD A, ONE ;CALCULATE 1-(X↑2)
03500 JUMPE A, ASIN1 ;WAS X EITHER -1.0 OR 1.0?
03600 PUSH P,A ;NO,
03700 PUSHJ P,$SQRT ;CALCULATE SQRT(1-X↑2)
03800 MOVE A,-1(P) ;GET THE ARGUMENT BACK AGAIN
03900 FDV A,B ;CALCULATE X/SQRT(1-X↑2)
04000 PUSH P,A ;THEN
04100 PUSHJ P,$ATAN ;CALCULATE ATAN(X/SQRT(1-X↑2)),.
04200 EXIT: SUB P,[XWD 2,2]
04300 JRST @2(P)
04400
04500 ASIN1: MOVE B, PIOT ;ANSWER IS EITHER PI/2 OR-PI/2
04600 SKIPG -1(P) ;WAS ORIGINAL ARGUMENT POSITIVE?
04700 MOVNS B ;NO, GET -PI/2
04800 JRST EXIT
04900
05000 TOOLRG: OUTSTR [ASCIZ /ASIN: ARGUMENT MAGNITUDE GREATER THAN 1.0; 0 RETURNED/]
05100 SETZ B,0 ;SET ANSWER TO ZERO.
05200 JRST EXIT ;RETURN
05300 PIOT: 201622077325 ;PI/2
05400 ONE: 1.0
05500
05600 BEND $ASIN
00100 ;FLOATING POINT SINGLE PRECISION ARCCOSINE FUNCTION
00200
00300 ;ACOS(X) IS CALCULATED IN THE FOLLOWING MANNER:
00400 ; IF X GREATER THAN 0, ACOS(X)=ATAN((SQRT(1-X↑2))/X)
00500 ; IF X LESS THAN 0, ACOS(X)=PI + ATAN((SQRT(1-X↑2))/X)
00600 ; IF X = 0, ACOS(X)=PI/2
00700
00800 ;THE RANGE OF DEFINITION FOR ACOS IS -1.0 TO +1.0.
00900 ;ARGUMENTS OUTSIDE OF THIS RANGE WILL CAUSE AN ERROR MESSAGE
01000 ;TO BE TYPED AND WILL RETURN AN ANSWER OF ZERO.
01100
01200 ;THE CALLING SEQUENCE FOR ACOS IS:
01300
01400 ; PUSH 17,ARG
01500 ; PUSHJ 17,$ACOS
01600 ;THE RESULT IS RETURNED IN AC 1
01700
01800
01900 BEGIN $ACOS
02000
02100 INTERNAL ACOS
02200
02300 ACOS:
02400 $ACOS: ;ENTRY TO ACOS ROUTINE.
02500 MOVM 1,-1(P) ;GET /ARG./ IN AC 1.
02600 CAMLE 1,ONE ;IS MAGNITUDE OF ARG. GT 1.0?
02700 JRST TOOLRG ;YES, GO TO ERROR RETURN.
02800 JUMPE 1,ZERARG ;IF ARG=0, GO TO ZERARG.
02900 FMPR 1,1 ;X↑2 IN AC 1.
03000 JFCL ;SUPPRESS ERROR MESSAGE FROM OVTRAP, IF NECESSARY.
03100 ;ON UNDERFLOW THE SPECIAL
03200 ;INTERRUPT CODE WILL SET
03300 ;AC 1 TO 0
03400 MOVNS 1 ;-X↑2 IN AC 1.
03500 FAD 1,ONE ;1.0-X↑2 IN AC 1.
03600 PUSH P,1
03700 PUSHJ P,$SQRT ;CALC. $SQRT(1.0-X↑2).
03800 FDVR 1,-1(P) ;($SQRT(1.0-X↑2))/X IS IN AC 1.
03900 JFCL ;SUPPRESS ERROR MESSAGE FROM OVTRAP, IF NECESSARY.
04000 ;ON UNDERFLOW THE SPECIAL INTERRUPT
04100 ;CODE WILL SET AC 1 TO 0
04200 ;ON OVERFLOW AC 1 WILL BE SET
04300 ;TO LARGEST MAGNITUDE WITH
04400 ;CORRECT SIGN
04500 PUSH P,1
04600 PUSHJ P,$ATAN ;FIND $ATAN($SQRT(1.0-X↑2)/X).
04700 SKIPL -1(P) ;SKIP IF ORIGINAL ARG LESS THAN 0.
04800 JRST EXIT ;RETURN.
04900 FAD 1,PII ;ANSWER IS PI + ANSWER IN AC 1.
05000 EXIT: SUB P,[XWD 2,2]
05100 JRST @2(P)
05200 ZERARG: MOVE 1,PI2 ;ANSWER IS PI/2.
05300 JRST EXIT ;
05400 TOOLRG: OUTSTR [ASCIZ /ACOS: ARGUMENT MAGNITUDE GREATHER THAN 1.0; 0 RETURNED/];
05500 SETZ 1,0 ;SET ANSWER TO ZERO.
05600 JRST EXIT
05700 ONE: 1.0
05800 PI2: 201622077325
05900 PII: 202622077325
06000
06100
06200 BEND $ACOS
00100 ;---------------------------------------------------------
00200
00300 ;PSEUDO RANDOM NUMBER GENERATOR AND INITIALIZING ROUTINE
00400 ;METHOD SUGGESTED BY D. H. LEHMER
00500
00600
00700 ;CALLING SEQUENCE FOR FUNCTION RAN:
00800
00900 ;PUSH 17,ARG
01000 ;PUSHJ 17,$RAN
01100 ;IF ARG NEQ 0, ARG IS USED AS PREVIOUS RANDOM NO.
01200 ;(I.E. AS THE STARTING VALUE), OTHERWISE THE PREVIOUS
01300 ;VALUE (WHICH WAS STORED LOCALLY) IS USED.
01400 ;ANSWER IS RETURNED IN ACCUMULATOR 1 AS A SINGLE
01500 ;PRECISION FLOATING POINT NUMBER IN THE RANGE
01600 ;0 LESS THAN X LESS THAN 1
01700
01800
01900 BEGIN $RAN
02000
02100 A←13
02200 B←14
02300
02400 INTERNAL RAN
02500
02600 RAN:
02700 $RAN:
02800 MOVE A,-1(P) ;IF ARG = 0 THEN
02900 JUMPE A,R1 ;USE PREVIOUS RANDOM NO.
03000 TLZ A,760000 ;OTHERWISE MASK 5 BITS
03100 MOVEM A,XN ;AND STORE NEW NO.
03200 R1: MOVE A,K ;GET K [14**29(MOD2**31 -1)]
03300 MUL A,XN ;MULTIPLY WITH LAST RANDOM NUMBER
03400 ASHC A,4 ;SEPARATE RESULT IN TWO 31 BIT WORDS
03500 LSH B,-4
03600 ADD A,B ;ADD THEM TOGETHER
03700 TLZE A,760000 ;SKIP IF RESULT LESS THAN 31 BITS
03800 ADDI A,1
03900 MOVEM A,XN ;STORE NEW RN IN INTEGER MODE
04000 HLRZ 1,A ;CONVERT TO FP IN TWO STEPS IN
04100 FSC 1,216 ;ORDER TO LOOSE NO LOW ORDER
04200 HRLI A,0 ;BITS
04300 FSC A,174
04400 FAD 1,A
04500 SUB P,[XWD 2,2]
04600 JRST @2(P) ;EXIT
04700 K: =630360016 ;14**29(MOD 2**31 -1)
04800 XN: =524287 ;STARTING VALUE
04900
05000 BEND $RAN
00100 ;---------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION HYPERBOLIC TANGENT ROUTINE
00300
00400 ;THIS ROUTINE CALCULATES THE TANH BY THE FOLLOWING ALGORITHM:
00500 ;IF ABSF(X) LESS THAN .00034, THEN TANH(X) = X
00600 ;IF ABSF(X) GREATER THAN 12.0, THEN TANH(X) = 1.0*SIGN(X)
00700 ;IF 0.17 LESS THAN= X LESS THAN 12.0, THEN TANH IS CALCULATED AS
00800 ; TANH(X) = 1.0 - 2(1.0 + EXP(2*X))**-1
00900 ;IF .00034 LESS THAN= X LESS THAN 0.17, THEN TANH IS CALCULATED AS
01000 ;TANH(X) = F(A+F↑2(B+C(D+F↑2)**-1))**-1
01100 ;WHERE X = 4*LOG(E) (BASE 2)
01200
01300 ;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
01400 ; PUSH 17,ARG
01500 ; PUSHJ 17,$TANH
01600 ;THE ANSWER IS RETURNED IN ACCUMULATOR 1
01700
01800 BEGIN $TANH
01900
02000 A← 1
02100 B←13
02200
02300
02400 INTERNAL TANH
02500
02600 TANH:
02700 $TANH: ;ENTRY TO TANH ROUTINE
02800 MOVE A,-1(P) ;PICK UP THE ARGUMENT
02900 MOVM B, A ;GET ABSF(ARGUMENT)
03000 CAMGE B, T1 ;RETURN TANH(X)=X IF
03100 JRST EXIT ;ABSF(X) .LE. .00034
03200 CAMLE B, T2 ;RETURN TANH(X) = 1.0*SIGN(X) IF
03300 JRST TH5 ;ARGUMENT GREATER THAN 12.0
03400 CAMGE B, T3 ;USE RATIONAL APPROXIMATION IF
03500 JRST TH3 ;ARGUMENT IS LESS THAN 0.17
03600 FMPRI B,202400 ;GET 2*ARG.
03700 PUSH P,B ;CALCULATE EXP(2X)
03800 PUSHJ P,$EXP
03900 MOVSI B, (1.0) ;FORM 1.0
04000 FAD A, B ;1 + EXP(2X)
04100 FDVM B, A ;(1 + EXP(2X))**-1
04200 FMPRI A,202400 ;2*(1 + EXP(2X))**-1
04300 FSBRM B, A ;1 - 2*(1 + EXP(2X))**-1
04400 SKIPGE -1(P) ;SKIP AHEAD IF ARG WAS GREATER THAN= 0.
04500 MOVNS A ;OTHERWISE,NEGATE THE ANSWER.
04600 EXIT: SUB P,[XWD 2,2]
04700 JRST @2(P)
04800
04900 TH3: FMP A, T7 ;FORM 4*X*LOG(E) BASE 2
05000 MOVEM A, TM1 ;SAVE IT IN TM1
05100 FMP A, A ;SQUARE IT
05200 MOVEM A, TM2 ;SAVE IT
05300 FAD A, T4 ;FORM F↑2 + T4
05400 MOVE B, T5 ;GET T5 IN ACCUMULATOR B
05500 FDV B, A ;T5/(F↑2 + T4)
05600 FAD B, T6 ;T6 + T5/(F↑2 + T4)
05700 FMP B, TM2 ;MULTIPLY BY F↑2
05800 FAD B, T7 ;ADD T7 (4*LOG(E) BASE 2)
05900 MOVE A, TM1 ;GET F IN ACCUMULATOR A
06000 TH5: FDV A, B ;DIVIDE F BY PARTIAL SUM
06100 JRST EXIT ;EXIT
06200
06300 T1: 165544410070 ;0.00034
06400 T2: 204600000000 ;12.0
06500 T3: 176534121727 ;0.17
06600 T4: 211535527022 ;349.6699888
06700 T5: 204704333567 ;14.1384514018
06800 T6: 173433723376 ;0.01732867951
06900 T7: 203561250731 ;5.7707801636
07000
07100 TM1: 0
07200 TM2: 0
07300
07400
07500 BEND $TANH
00100 ;---------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION HYPERBOLIC COSINE FUNCTION.
00300
00400 ;COSH(X) IS CALCULATED AS FOLLOWS:
00500 ; IF ABS(X) LESS THAN= 88.029,
00600 ; COSH(X) = 1/2(EXP(X) + 1.0/EXP(X))
00700 ; IF ABS(X) GREATER THAN 88.029 AND (ABS(X)-LN(2)) LESS THAN= 88.029,
00800 ; COSH(X) = EXP(ABS(X)-LN(2))
00900 ; IF (ABS(X)-LN(2)) GREATER THAN 88.029,
01000 ; COSH(X)=377777777777
01100 ; AND AN ERROR MESSAGE IS RETURNED.
01200
01300 ;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
01400 ; PUSH 17,ARG
01500 ; PUSHJ 17,$COSH
01600 ;THE ANSWER IS RETURNED IN AC 1.
01700
01800 BEGIN $COSH
01900
02000 A←13
02100
02200
02300 INTERNAL COSH
02400
02500 COSH:
02600 $COSH: ;ENTRY TO HYPERBOLIC COSINE ROUTINE.
02700 MOVE 1,-1(P) ;PICK UP THE ARGUMENT.
02800 MOVM A,1 ;PUT ABS(X) IN A
02900 CAMLE 2,EIGHT8 ;IF ABS(X) GREATER THAN 88.029,
03000 JRST OV88 ;GO TO OV88.
03100 PUSH P,A ;O'E, CALC. EXP(ABS(X))
03200 PUSHJ P,$EXP
03300 MOVSI A,(1.0) ;PUT 1.0 IN A
03400 FDVR A,1 ;CALC. 1.0/EXP(ABS(X)).
03500 FADR 1,A ;CALC. EXP(ABS(X)) + EXP(-ABS(X)).
03600 FDVRI 1,202400 ;DIVIDE THIS BY 2.0.
03700 EXIT: SUB P,[XWD 2,2] ;RETURN.
03800 JRST @2(P)
03900
04000 OV88: FSBR A,LN2BE ;FORM ABS(X)-LN(2).
04100 CAMG A,EIGHT8 ;OVERFLOW?
04200 JRST EXPP ;NO,GO AHEAD.
04300 OUTSTR [ASCIZ /COSH: RESULT TOO LARGE - LARGEST POSITIVE NUMBER RETURNED./];
04400 HRLOI 1,377777 ;ANSWER = +INFINITY.
04500 JRST EXIT ;RETURN
04600
04700 EXPP: PUSH P,A ;CALC. EXP(ABS(X)-LN(2)).
04800 PUSHJ P,$EXP
04900 JRST EXIT ;RETURN.
05000
05100 EIGHT8: 207540074636 ;88.029
05200 LN2BE: 200542710300 ;LOG(2) BASE E.
05300
05400 BEND $COSH
00100 ;---------------------------------------------------------
00200 ;FLOATING POINT SINGLE PRECISION HYPERBOLIC SINE FUNCTION.
00300
00400 ;SINH IS CALCULATED AS FOLLOWS:
00500 ; IF ABS(X) GREATER THAN 88.029,
00600 ; SINH(X)=(EXP[ABS(X)-LN(2)])*SIGN(X)
00700 ; IF ABS(X) LESS THAN= 0.10,
00800 ; SINH(X)=X+(X**3)/6+(X**5)/120
00900 ; FOR ALL OTHER VALUES OF X,
01000 ; SINH(X)=1/2[EXP(X)-1/EXP(X)]
01100
01200 ;THE CALLING SEQUENCE IS:
01300 ; PUSH 17,ARG
01400 ; PUSHJ 17,$SINH
01500
01600 ;THE ANSWER IS RETURNED IN AC 1.
01700
01800 BEGIN $SINH
01900
02000 A←13
02100 B←14
02200
02300 INTERNAL SINH
02400
02500 SINH:
02600 $SINH: ;ENTRY TO HYPERBOLIC SINE ROUTINE.
02700 MOVE 1,-1(P) ;PICK UP THE ARG.
02800
02900 MOVM A,1 ;GET MAGNITUDE OF ARG IN A
03000 CAMLE A,EIGHT8 ;IF ABS(X) GREATER THAN 88.029,
03100 JRST OV88 ;THEN GO TO OV88.
03200 CAMG A,ONE10T ;IF ABS(X) LESS THAN= 0.10,
03300 JRST SERIES ;THEN GO TO SERIES.
03400 PUSH P,A ;CALCULATE EXP(ABS(X)).
03500 PUSHJ P,$EXP ;ABS(X) IS IN A
03600 HRLZI B,576400 ;PUT -1.0 IN B
03700 FDVR B,1 ;CALC. -EXP(-ABS(X)).
03800 FADR 1,B ;CALC. EXP(ABS(X))-EXP(-ABS(X)).
03900 FDVRI 1,202400 ;CALC. THIS/2.0
04000 SKIPGE -1(P) ;ANSWER IS POSITIVE.
04100 MOVNS 1,1 ;ANSWER IS NEGATIVE.
04200 EXIT:
04300 SUB 17,[XWD 2,2]
04400 JRST @2(P)
04500 SERIES: FMPR A,A ;CALC. X↑2.
04600 JFCL ;SUPPRESS ERROR MESSAGE FROM OVTRAP.
04700 ;ON UNDERFLOW, SPECIAL INTERRUPT
04800 ;CODE RETURNS 0.
04900 MOVEM A,SX2 ;SAVE X↑2 IN SX2.
05000 FDVR 2,ONE120 ;CALC.X↑2/120
05100 JFCL ;SUPPRESS ERROR MESSAGE FROM OVTRAP.
05200 ;ON UNDERFLOW, SPECIAL INTERRUPT
05300 ;CODE RETURNS 0.
05400 FADR A,ONESIX ;CALC. (X↑2/120)+1/6
05500 FMPR A,SX2 ;MULTIPLY IT BY X↑2.
05600 JFCL ;SUPPRESS ERROR MESSAGE FROM OVTRAP.
05700 ;ON UNDERFLOW SPECIAL INTERRUPT
05800 ;CODE RETURNS 0.
05900 FADRI A,(1.0) ;ADD 1.0.
06000 FMPR 1,A ;MULTIPLY BY X.
06100 JRST EXIT ;RETURN.
06200 OV88: FSBR A,LN2BE ;CALC.ABS(X)-LN(2)
06300 CAMG A,EIGHT8 ;OVERFLOW?
06400 JRST EXPP ;NO,GO TO CALC.
06500 OUTSTR [ASCIZ /SINH: RESULT TOO LARGE - LARGEST POSITIVE NUMBER RETURNED/];
06600 HRLOI 1,377777 ;SET ANS.=INFINITY.
06700 JRST EXPP+2 ;GO TO SET SIGN OF ANS.
06800
06900 EXPP: PUSH P,A ;CALC. EXP
07000 PUSHJ P,$EXP
07100 SKIPGE -1(P) ;RETURN ANS. GREATER THAN 0 IF X GREATER THAN 0.
07200 MOVNS 1,1 ;O'E, ANS. LESS THAN 0.
07300 JRST EXIT ;RETURN.
07400
07500 LN2BE: 200542710300 ;LN(2)
07600 EIGHT8: 207540074636 ;88.029
07700 ONE10T: 0.10
07800 SX2: 0
07900 ONE120: 207740000000 ;120.0
08000 ONESIX: 0.16666667
08100
08200 BEND $SINH
00100 ;---------------------------------------------------------
00200
00300 ;FLOATING POINT SINGLE PRECISION ARCTANGENT OF TWO ARGUMENTS
00400 ;RETURNS ARCTANGENT OF A/B
00500 ;IF ARGUMENT IS IN 2ND QUADRANT, ATAN2(A/B) = PI + ATAN(A/B)
00600 ;IF ARGUMENT IS IN 3RD QUADRANT, ATAN2(A/B) = ATAN(A/B) - PI
00700 ;IF A/B OVERFLOWS (OR DIVIDE CHECK), THEN RETURN
00800 ; +PI/2 IF A GREATER THAN= 0, AND
00900 ; -PI/2 IF A LESS THAN 0.
01000 ;IF A/B UNDERFLOWS, THEN RETURN
01100 ; 0 IF B GREATER THAN= 0, AND
01200 ; +PI IF B LESS THAN 0 AND A GREATER THAN= 0,
01300 ; -PI IF B LESS THAN 0 AND A LESS THAN 0.
01400
01500 ;THERE IS NO RESTRICTION ON THE ARGUMENTS
01600
01700 ;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:
01800 ; PUSH 17,ARG1
01900 ; PUSH 17,ARG2
02000 PUSHJ 17,$ATAN2
02100 ;THE ANSWER IS RETURNED IN ACCUMULATOR 1.
02200
02300 BEGIN $ATAN2
02400
02500 A←13
02600 B←1
02700
02800 INTERNAL ATAN2
02900
03000 ATAN2:
03100 ↑$ATAN2: ;ENTRY POINT TO ATAN2 ROUTINE
03200 MOVE A,-2(P) ;PICK UP FIRST ARGUMENT
03300 MOVE B,-1(P) ;PICK UP SECOND ARGUMENT
03400 FDVR A, B ;FORM A/B
03500 JFCL 0,OVUNFO ;SUPPRESS ERROR MESSAGE FROM
03600 ;OVTRAP IF NECESSARY AND GO TO
03700 ;OVUNFO IF AN EXCEPTION OCCURRED
03800 ;SPECIAL INTERRUPT CODE SETS
03900 ;OVPCWD AND RETURNS BEST VALUE
04000 ;IT CAN IN A.
04100 PUSH P,A ;CALCULATE ATAN(A/B)
04200 PUSHJ P,$ATAN
04300 SKIPL -1(P) ;IF B GREATER THAN 0, SGN(ATAN2)=SGN(A)
04400 JRST EXIT ;EXIT
04500 JUMPGE B, ATAN2A ;IS B POSITIVE?
04600 FADR B, PII ;NO, SECOND QUADRANT, ADD PI
04700 EXIT: SUB P,[XWD 3,3]
04800 JRST @3(P)
04900
05000 ATAN2A: FSBR B, PII ;YES,3RD QUADRANT, SUBTRACT PI
05100 JRST EXIT ;EXIT
05200 OVUNFO: MOVE A,OVPCWD ;PICK UP FLAGS.
05300 TLNE A,000100 ;SKIP IF OVERFLOW.
05400 JRST UNDER ;O'E GO TO UNDERFLOW CASE.
05500 MOVE B,HALFPI ;ANSWER SET TO PI OVER TWO.
05600 SKIPGE -2(P) ;SKIP IF ANS IS TO BE POSITIVE.
05700 MOVNS B ;SET ANSWER NEGATIVE.
05800 JRST EXIT ;EXIT.
05900
06000 UNDER: JUMPL B,BNEG ;IF B GREATER THAN= 0
06100 MOVEI B,0 ;THEN RETURN 0
06200 JRST EXIT
06300 BNEG: MOVE B,PII ;ELSE IF B LESS THAN 0
06400 SKIPGE -2(P) ;RETURN PI IF A GREATER THAN= 0
06500 MOVNS B ;OR -PI IF A LESS THAN 0
06600 JRST EXIT
06700
06800 PII: 202622077325 ;3.141592653589793238462643...
06900 HALFPI: 201622077325
07000
07100 BEND $ATAN2
07200
07300
07400 XLIST ;NO SEE LITERALS
07500 END
07600
07700